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Abstract 



We present a number of complete sets of series expansion formulas for neutrino oscillation 
probabilities in matter of constant density for three flavors. In particular, we study expansions 
in the mass hierarchy parameter a = Am'^i/Am'^^ and mixing parameter S13 = sin6'i3 up to 
second order and expansions only in a and only in S13 up to flrst order. For each type of 
expansion we also present the corresponding formulas for neutrino oscillations in vacuum. We 
perform a detailed analysis of the accuracy of the different sets of series expansion formulas 
and investigate which type of expansion is most accurate in different regions of the parameter 
space spanned by the neutrino energy E, the baseline length L, and the expansion parameters 
a and S13. We also present the formulas for series expansions in a and in S13 up to flrst order 
for the case of arbitrary matter density proflles. Furthermore, it is shown that in general all 
the 18 neutrino and antineutrino oscillation probabilities can be expressed through just two 
independent probabilities. 
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1 Introduction 



The discovery of neutrino oscillations in atmospheric, solar, and reactor neutrino experiments 
has turned neutrino physics into one of the most exciting and active fields of particle physics. 
By now, a significant amount of information on neutrino properties has been obtained. The 
results of the atmospheric neutrino experiments [1-3] and the K2K accelerator neutrino ex- 
periment [4] have allowed the determination of the fundamental neutrino parameters Amgj^ 
and ^23 to an accuracy of about 30% and 15%, respectively, while the solar neutrino experi- 
ments [5] and the KamLAND reactor neutrino experiment [6] have measured the parameters 
Am|]^ and 612 to an accuracy of about 15%. In fact, this means that neutrino physics is 
now entering an era of precision measurements of the neutrino oscillation parameters. Future 
experiments with superbeams and neutrino factories will determine the "atmospheric" and 
"solar" neutrino oscillation parameters to an accuracy of the order of 1%. They are also 
expected to measure the elusive leptonic mixing angle 6*13 for which at present only an upper 
limit exists [7] , or to put a more stringent limit on this angle, as well as to clarify the issues 
of the neutrino mass hierarchy and possibly of leptonic CP violation. Important information 
on neutrino oscillation parameters can also be obtained from the future atmospheric, solar, 
reactor, and supernova neutrino experiments. 

Increasing accuracy and reach of the present and especially of forthcoming experiments 
put forward new and more challenging demands to the theoretical description of neutrino 
oscillations. In order to be able to determine the fundamental neutrino parameters from the 
data with high precision, one needs, among other things, very accurate theoretical expressions 
for the probabilities of neutrino oscillations in matter and in vacuum. While these proba- 
bilities can, in principle, be calculated numerically with any requisite accuracy, it is highly 
desirable to have also analytic expressions for them. Such analytic expressions would reveal 
the basic dependence of the neutrino oscillation probabilities on the fundamental neutrino 
parameters and on the characteristics of the experiment and thus facilitate the choice of the 
experimental setup as well as the analysis of the data. They would also help to understand 
the physics underlying various flavor transitions and to resolve the parameter degeneracies 
and other ambiguities, such as fundamental versus matter-induced CP violation. 

The purpose of this paper is to present a collection of approximate analytic formulas 
for the neutrino oscillation probabilities. We derive a number of complete sets of series 
expansion formulas for the three-flavor neutrino oscillation probabilities in matter and in 
vacuum. The probabilities are expanded in the mass hierarchy parameter a = Aml^/Aml^, 
the mixing parameter S13 = sin^^i3, or in both of them. We also study in detail the accuracy 
of the obtained expressions in different regions of the parameter space and identify the "best 
choice" in each case of interest. 

Before proceeding to present our results, we give here a brief overview of the previous 
work on the subject. Analytic formulas for three-flavor neutrino oscillation probabilities have 
been derived in a number of papers. Exact formulas for the neutrino oscillation probabilities 
in vacuum can be found, e.g., in Ref. [8]. Exact expressions can also be obtained in the 
case of three-flavor neutrino oscillations in matter of constant density [9-16]. However, the 
corresponding formulas are rather complicated and not easily tractable. This also applies 
to the exact analytic three-flavor formulas obtained for some special cases of non-uniform 
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matter density: linear matter density [17] and exponentially varying matter density [18]. 

Approximate analytic formulas for three-fiavor neutrino oscillation probabilities in mat- 
ter have been derived in a number of papers, see e.g., Refs. [19-24]. In Refs. [19,20], the 
Mikheyev-Smirnov-Wolfenstein (MSW) resonances [25] in matter of varying density were 
studied assuming large separation and independence of the high-density and low-density res- 
onances. In Refs. [20,21,23], the adiabatic approximation was used to derive three-flavor neu- 
trino oscillation probabilities in matter of varying density, whereas in Rcf. [24], the Magnus 
expansion for the time evolution operator was used, which is equivalent to the anti-adiabatic 
approximation. 

In a number of papers an approach similar to ours was adopted, i.e., the neutrino os- 
cillation probabilities were expanded in the small parameters a, Sis, or in both of them. 
In Ref. [26], exact analytic expressions for three-flavor neutrino oscillation probabihties in 
matter with an arbitrary density profile were obtained in the limit a — >^ by reducing the 
problem to an effective two-flavor one. In Ref. [27], a similar approach was employed to ob- 
tain the neutrino oscillation probabilities in matter of arbitrarily varying density in the limit 
^13 0, whereas in Refs. [28,29], expressions up to first order in S13 were derived. For the 
case of matter of constant density, the hmit a — > was considered in Ref. [30]. Expansions 
up to first order in a were carried out in Refs. [31.32]. In Refs. [33,34], both the solar mass 
squared difference lS.vn^^ and matter effects were treated as perturbations and the transition 
probabilities up to first order in them were derived. Expansions in both a and S13 were used 
in Refs. [35-38] . We note that the neutrino oscillation probabilities derived in most of the 
above-mentioned papers either did not constitute a complete set [i.e., a set of probabilities 
from which the probabilities in all channels can be obtained), or contained expressions from 
which some terms were missing, especially in the case of the probabilities Pivpu t'^) and 
Piy\x ^t)- Thus, to the best of our knowledge, our study is the first one in which, for the 
case of matter of constant density, complete and consistent expansions in a and 613 to second 
order and expansions only in a and only in S13 up to first order are performed. 

The paper is organized as follows. In Sec. 2 we set the general formalism and notation, 
discuss some general relations satisfied by the oscillation probabilities, and show that all 
the 18 neutrino and antineutrino oscillation probabilities can be expressed through just two 
independent probabilities. In Sec. 3 we present a series expansion of the neutrino oscillation 
probabilities up to second order in both a and S13 (the so-called "double expansion"), whereas 
in Sees. 4 and 5 we consider the probabilities expanded up to first order in a and S13, 
respectively ("single expansions"). In all three cases, we give the probabilities for matter of 
constant density (Sees. 3.1, 4.1, and 5.1) and in vacuum (Sees. 3.2, 4.2, and 5.2). We also 
compare our formulas with the corresponding expressions existing in the literature (when 
available), pointing out agreements and disagreements. In Sec. 6, we discuss the qualitative 
behavior of the neutrino oscillation probabilities, the relevance of matter effects, and give a 
detailed evaluation of the accuracy of the various formulas. Furthermore, we comment on the 
application of our probability formulas to neutrino oscillation experiments. We summarize 
our results in Sec. 7. Finally, several methods that have been used to derive the formulas are 
presented in the appendices. In Appendix A, we describe the perturbative diagonalization 
of the effective Hamiltonian of the neutrino system in the case of matter of constant density, 
while in Appendix B, the details of the perturbative expansion of the neutrino evolution 
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equation in the case of arbitrary matter density profiles are given. 



2 General formalism and notation 

We consider three-flavor neutrino oscillations and adopt the standard parameterization of the 
leptonic mixing matrix U [39]: 

U = 023UsOisUlOu 

(C12C13 S12C13 si3e~"''^p\ 

-^12^23 - Ci2Si3S23C"^^P C12C23 - Sl2Sl3S23C"^'^'' C13S23 • (l) 

S12S23 - Ci2Sl3C23e"^'^'' -C12S23 - Sl2Sl3C23e"^'^P C13C23 / 

Here Oij is the orthogonal rotation matrix in the ij'-plane which depends on the mixing 
angle 6ij, Us = diag(l, 1, e"^c;p), 6cp being the Dirac-type CP-violating phase, Sij = sinOij 
and Cij = cos 9ij. In the three-flavor case there are also, in general, two Majorana-type CP- 
violating phases; however, these phases do not affect neutrino oscillations, and therefore will 
not be considered here. Without loss of generality, one can assume all the mixing angles to 
lie in the first quadrant {i.e., between and 7r/2), while the CP-violating phase 5cp is allowed 
to lie in the interval [0, 2n]. 

Let us denote by Pap = Pi^a ~^ '^/s) the probability of transition from a neutrino flavor 
a to a neutrino flavor f3, and similarly for antincutrino flavors, i.e., P^p = P{ua ^13) ■ 
In general, the three-flavor neutrino oscillation probabilities in matter depend on eight 
parameters and one function: 

Pa3 = Pa(3{^ml^,^ml^,0i2, 9i3,923, Sep; E,L,V{x)), a,(3 = e,(i,T. (2) 

Here Am|^- = mf — th?- are the neutrino mass squared differences, E is the neutrino en- 
ergy, L is the baseline length, and V{x) is the matter-induced effective potential, x G [0, L] 
being the coordinate along the neutrino path. The neutrino mass squared differences, the 
leptonic mixing angles, and the CP-violating phase are fundamental parameters and thus 
experiment-independent, whereas the neutrino energy, the baseline length, and the matter- 
induced effective potential vary from experiment to experiment. The present best-flt values 
and 3(7 allowed ranges of the fundamental neutrino parameters found in a recent global flt 
of the neutrino oscillation data [40] are summarized in Table 1. Unless otherwise stated, 
all calculations in the present paper are performed for the following values of the neutrino 
parameters: ^m^i = 7 ■ 10~^eV^, 612 — 33°, and 623 = 45°. For the atmospheric mass 
squared difference Amg^ we adopt the current best-flt value given by the Super-Kamiokande 
Collaboration, |Am|^| = 2 ■ 10^^ eV^ [2], which is slightly smaller than the value given in Ta- 
ble 1. The sign of Am|^ is related to the neutrino mass hierarchy: for the normal (inverted) 
hierarchy one has Amg^ > (Am|i < 0). For the leptonic mixing angle ^13 we allow values 
below the 90% C.L. upper bound found in the global fit of the neutrino oscillation data [40] 
for [Am^^l fixed at 2 • lO'^eV^: 

ei3<10.8°, or si3 = sin^i3 <0.19, or sin^ 2^13 < 0.14 . (3) 

For the CP violation phase 6cp we allow values between and 2tt. 
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Parameter 


Best-fit value 


Range [6a) 




6.9 ■ 10-5 cV^ 


(5.4^9.5) ■ 10-5 cV^ 




2.6 ■ 10-3 eV^ 


(1.5^3.7) ■ 10-3 eV^ 




33.2° 


28.6° 38.6° 




4.4° 


^ 13.4° 


023 


46.1° 


33.8° ^58.1° 


Sep 




0^27r 



Table 1: Present best-fit values and 3a allowed ranges of the fundamental neutrino param- 
eters from a three-flavor flt to global neutrino oscillation data [40] . 



Inspecting the values of the fundamental neutrino parameters in Table 1, one can identify 
two natural candidates for small expansion parameters of the neutrino oscillation probabil- 
ities. These are the small leptonic mixing angle ^13 (or, equivalently, S13) and the mass 
hierarchy parameter 



a = -— |i ~ 0.026 , (0.018) 0.021 <a< 0.036 (0.053) at 90% C.L. {3a) , (4) 

where we have taken the current best-fit value and allowed ranges from Ref. [40]. In the 
following, we will derive a number of formulas for series expansions of the neutrino oscillation 
probabilities in these small quantities. Comparing Eqs. (3) and (4), one realizes that current 
data constrain the parameter a to a relatively narrow range, while S13 is only bounded from 
above and might be significantly larger or smaller than a. The relative size of these two 
expansion parameters will be important for the validity of a given type of expansion. 

In order to find the neutrino oscillation probabilities for a given experimental setup, one 
has, in general, to solve the Schrodinger equation for the neutrino vector of state in the fiavor 
basis |z/(t)) = {Ue{t) u^{t) i'r{t)Y'- 

■4\v{t)) = H\v{t)) (5) 



with the effective Hamiltonian 



H ~ ^t/diag(0, Ami,, ^rnl,)U^ + diag(V, 0, 0) . (6) 

Here V is the charged-current contribution to the matter-induced effective potential of v^, [25] . 
We have disregarded the neutral-current contributions to the neutrino potentials in matter, 
since they are the same for z/g, z/^, and ^ and so do not affect neutrino oscillations. Note that 
Eq. (6) holds for neutrinos, whereas for antineutrinos one has to perform the replacements 

U ^U* , V ^-V . (7) 



^Up to tiny radiative corrections [41] which are neghgible except at very high densities available, e.g., 
inside supernovae. 
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The potential V{x) is given in convenient units by 

V(x) ~ 7.56 X 10"^^ ( -4^^ Yeix) eV , (8) 
Vg/cmV 

where p{x) is the matter density along the neutrino path and Ye{x) is the number of electrons 
per nucleon. For the matter of the Earth one has, to a very good accuracy, Ye ~ 0.5. 

For many practical applications (such as long-baseline accelerator experiments, as well 
as oscillations of atmospheric, solar, and supernova neutrinos inside the Earth when they 
do not cross the Earth's core) it is a very good approximation to assume that the matter 
density along the neutrino trajectory is constant (see, e.g., Refs. [42-44]). Typical values for 
the matter density are Pcmst — 3g/cm^ in the Earth's crust and Pmantie — 4.5g/cm^ in its 
mantle. In situations where the neutrinos also cross the Earth's core or for strongly varying 
matter density profiles like those inside the Sun or supernovae, the constant matter density 
approximation is not valid. 

The neutrino oscillation probabilities can be found as Pap — \Si3a{t,tQ)\^, where S{t,to) 
is the evolution matrix such that 

\iy{t)) = S{t,toMto)) , S{to,to) = l. (9) 

Note that S{t,tQ) satisfies the same Schrodinger equation, Eq. (5), as \iy{t)). In the case 
of matter of constant density, the evolution matrix can be obtained by diagonalizing the 
Hamiltonian in Eq. (6) according to H — U'HU'\ where U' is the leptonic mixing matrix in 
matter and H — diag(£'i, £'2, E^). The evolution matrix is then given by 

3 

^/3a(i, io) = Y.^Kd*U'p,e-''^^'^ , a, /3 = e, T , (10) 

i=l 

where we have identified L = t — tg. 

Before presenting our results for the neutrino oscillation probabilities, we discuss some of 
their general properties as well as relations between them. First, we note that Eq. (7) implies 
that one can relate the oscillation probabilities for antineutrinos to those for neutrinos by 

Pais = PapiScp ^ -Sep, V^-V), a,f3^e,fi,T. (11) 

Second, in general, i.e., both in vacuum and in matter with an arbitrary density profile, it 
follows from the unitarity of S{t, to) (conservation of probability) that 

^Pa/s = '^Pap ^i, a,(3^e,ii,T. (12) 

a /? 

These relations imply that five out of the nine neutrino oscillation probabilities can be ex- 
pressed in terms of the other four [45] . 

Besides these general properties, there exists an additional symmetry due to the specific 
parameterization of the mixing matrix given in Eq. (1) and the fact that the rotation matrix 
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O23 commutes with the matter potential term of the Hamiltonian in Eq. (6). It is easy to 
show that the evolution matrix can be written as 



S{t,to)^02sS'{t,to)Ol (13) 

where S'{t,to) does not depend on 623. This can be used to prove some useful relations 
between the probabihties. Let us denote ^ 

= ^0/3(^23 ^ C23 , sin 26^23 -sin 26^23) , a,P = e,iJ,,T. (14) 
Using Eqs. (13) and (14), one can readily show that 

while Pee turns out to be independent of ^23- 

Out of the three conditions in Eq. (15), only two are independent, as each of them 
can be derived from the other two and the unitarity conditions (12). Hence, the number 
of independent neutrino oscillation probabilities is reduced to two. Thus, we come to the 
important conclusion that all the nine neutrino oscillation probabilities can be expressed 
through just two independent probabilities provided that their dependence on the mixing 
angle ^23 is known. However, the choice of these independent probabilities is restricted: 
they should not include Pee; which is independent of ^23; nor should they be a pair of the 
probabilities, which are T-reverse of each other, or go into each other (or T-reverse of each 
other) under the transformation S23 ^ C23 > sin 2^23 — — sin 2^23- 

One possible choice, and the one that we will use, is Pe/j, and P^,-. For completeness, we 
give here the expressions for all the other neutrino oscillation probabilities in terms of these 
two. Using Eqs. (12) and (15), one easily finds 

Pee = l-(Pe;. + Pe;.), (16) 

Per Pefi ; (1''') 

Pfj,e PeiJ, P/XT ~l~ Pfj,T -I (-i^) 

Pjj-ix 1 Pefx PfjiT 1 (-^'^) 

Pre Pefi ^~ P^t P^t i (20) 

Pth — PfXT -I (21) 

P,, = l-(Pe^ + P^,). (22) 

In addition to the above relations, one can study the transformations of the neutrino 
oscillation probabilities under the time reversal Pap —>■ Ppa- It can be shown [28] that in 
matter with an arbitrary density profile 



13a 



PapiScp -Sep, V{x) Vrey{x)) , a,(3^e,ii,T, (23) 



^The transformation in Eq. (14) can be achieved, e.g., through the substitution 623 — > O23 + 7r/2 or 
6*23 ^ 6*23 + 37r/2. 
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where Vrcv(3^) is the "reverse" potential, which corresponds to the interchanged positions 
of the neutrino source and the detector. In the case of symmetric matter density profiles 
(including matter of constant density), Kev(a^) = y{x), and Eq. (23) simplifies to 

Pf}a ^ Pa(}{Scp ^ -Sep) : a,p^e,iJ,T. (24) 

While Eq. (24) does not further reduce the number of independent probabihties, it yields 
relations which can be useful for cross-checking the formulas for Pq,^ in the case of matter 
with symmetric density profiles. 

By applying the rule given in Eq. (11), one can obtain from Eqs. (16)-(22) the cor- 
responding probabilities for the antineutrino oscillations. Thus, the expressions for all 18 
probabilities of neutrino and antineutrino oscillations can be found from the formulas for 
just two independent neutrino oscillation probabilities, which, as was already mentioned, we 
choose to be Pg^ and P^^. In order to be more explicit, we will in some cases also give 
formulas for additional neutrino oscillation channels. 

In the following Sees. 3, 4, and 5, wc give our results for the various series expansions of 
the neutrino oscillation probabilities. We will adopt the following abbreviations: 

Am|i 2A ^ ' 



3 Series expansion up to second order in a and sis 



3.1 Matter of constant density 



In this section, we present the series expansion formulas for three-flavor neutrino oscilla- 
tion probabilities in matter of constant density up to second order in both a and S13. The 
probabilities are calculated by diagonalizing the Hamiltonian (6) up to second order in these 
parameters, as described in Appendix A.l. We find for the eigenvalues of the Hamiltonian 

Ami-, , o A ,sin^2^i2\ , , 

S,.^(«4-.= !5^). (28) 

Am?! / ^ A \ , , 

Calculating the mixing matrix in matter U' as described in Appendix A.l and using Eq. (10) 
for the evolution matrix S, it is straightforward (although somewhat tedious) to derive the 
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following expressions for the neutrino oscillation probabilities: 

n n ^ sin^AA , sin^(A-l)A , , 

Pee = 1 - a' sin^ 2^12 - 4 , (30) 

2.2^^ 2 sin^^A ^ 2 2 sin^(A-l)A 

= Sm^ 2^12 C^3 ^2 + ^ 4 _ ^^2 

. . /AT , sinAA sin(A - 1)A 
+ 2 a si3 sm 2^12 sm 2^23 cos(A - 5cv) A- I ' ^ ^ 

2-2^/1 2 sinMA ^ 2 2sin^(^-l)^ 

Per = Sm^ 2^12 ^23 ^2 + ^ S?3 4 _ -^^2 

, sinAA sin(A - 1)A , , 

- 2 q; si3 sin 2^i2 sin 2^23 cos(A - 5cp) A- I ' ^^^^ 

P^^ = 1 - sin^ 26*23 sin^ A + o; sin^ 26*23 A sin 2 A 

sin^ AA 

- sin^ 26*12 C23 — C12 sin^ 26*23 cos 2 A 

1 9 9 ^ 9 ^ / . sin A A , , , . A 
+ — sin^ 26*12 sin^ 26*23 sin A — - — cos(^ - 1) A - — sin 2A 

_ sin^(A-l)A 

4^13 523 (^_^)2 

2 9 •) n f . . . sin(.4 - 1)A A , 

sm^ 26*23 sm A cos AA — A sm 2A 



A-l ^ - 1 2 

sin^lA sin(A- 1)A 



2 a S13 sin 26*12 sin 26*23 cos 5cp cos A 



A A-l 



+ — - q;si3 sin 2^12 sin 26*23 cos 26*23 cos^cp sin A ^74sin A — ^—^ — cos(74 — 1)A^ , 

(33) 

P^r = sin^ 26*23 sin^ A - a C12 sin^ 26*23 A sin 2 A + sin^ 26*23 A^ cos 2 A 

1 9.9^^ .9^^ f . A sin A A , , , , . A . ^ . 
- TttoT sm^ 26*12 sm^ 26*23 sm A — - — cos(^ - 1)A - — sm2A 

2 9 . 9^^ / . A .Asin(A-l)A A^ . 
+ sm^ 26*23 sm A cos AA — ^- —A sm 2 A 

A sinAA sin(A - 1)A 
+ 2 a Si3 sm 26*i2 sm 26*23 sm ocp sm A — 
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2 / sin ^4 A \ 
— — - a Si3 sin 2^12 sin 2^23 cos 2^23 cos 5cp sin A (A sin A — — cos{A — 1) A j . 

(34) 



Formally, our calculations are based upon the approximations a, S13 -C 1 and no explicit 
assumptions about the values of L/E are made. However, we remark that the series expansion 
formulas (30)-(34) are no longer valid as soon as a A = AmliL/^iE) becomes of order unity, 
i.e., when the oscillatory behavior due to the "solar" mass squared difference Am|i becomes 
relevant. This can happen for very long baselines and/or very low energies. See also Sec. 6 
for a detailed discussion of the accuracy of these formulas. 

Prom Eqs. (27)-(29) one can see that in vacuum {A — 0) and at the atmospheric reso- 
nance {A = 1) the expressions for the eigenvalues are divergent, and one would expect the 
expansion to break down. In these cases, two out of the three eigenvalues of the unperturbed 
Hamiltonian are degenerate and, strictly speaking, the degenerate perturbation theory rather 
than the ordinary one should be employed. However, it turns out that, though the eigenvalues 
(27)-(29) are divergent, the neutrino oscillation probabilities are finite in the limits A — > 
and A — >■ 1. The reason for this interesting behavior is a cancellation of divergences between 
the eigenvalues and the matrix elements of the leptonic mixing matrix in the calculation of 
the evolution matrix according to Eq. (10). In particular, in the limit /I — > 0, the correct 
vacuum neutrino oscillation probabilities are obtained. 

We shall now compare the above results with those existing in the literature. Equa- 
tions (27)-(32) have previously been derived in Ref . [35] and confirmed in Ref . [36] . Expres- 
sions (33) and (34) arc new. In Ref. [38], an expression for P^^. was found to first order in 
a, including the 0{a S13) term, which can be compared with the corresponding terms in our 
Eq. (34). We find that, while our C(l) and 0{a) terms coincide with those in Eq. (A3) of 
Ref. [38], our 0{asi3) term is quite different. In particular, we disagree with the statement 
in Ref. [38] that to order a the probabihty Pfj,r does not depend on the CP- violating phase 
5cp- The existence of a term proportional to asin^cp in P^^- is actually expected, since in 
matter of constant density the phase 6cp is the sole source of T-violation, and up to the sign, 
the T-odd terms in all three transition probabilities must be the same [26]. This is indeed 
seen in Eqs. (31), (32), and (34). We also note that the term of order asia in Eq. (A3) of 
Ref. [38] diverges at the atmospheric resonance (A — l), while our expression (34) is regular 
at all physical values of parameters. 

3.2 Vacuum neutrino oscillation probabilities up to second order in oc and S13 

The vacuum neutrino oscillation probabilities up to second order in a and S13 are given by 

p;^ = 1 - q;2 sin2 2^12 A^ - 4 siu^ A , (35) 

P^^ = cc^ sin^ 26112 C23 + A sj^ sj^ sin^ A + 2asi3 sin 26I12 sin 26'23 cos( A - 5cp) A sin A , 

(36) 

P^ = sin^ 2^12 S23 A^ + 4 5^3 C23 sin^ A - 2 a S13 sin 2^12 sin 2^23 cos (A - (5cp) A sin A , 

(37) 
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P^^Jf = 1 - sin^ 2^23 sin^ A + a c?2 sin^ 2^23 A sin 2 A 



MM 



- A^ [sin^ 26I12 C23 + cl^ sin^ 26I23 (cos 2 A - S12)] + 4 s% S23 cos 26'23 sin^ A 

— 2 q; Si3 sin 2^12 S23 sin 2^23 cos 5cp A sin 2 A , 

P,!r = sin^ 2^23 sin^ A - a c^2 sin^ 2^23 A sin 2 A 



+ a' sin^ 2^23 A^ \ c\^ cos 2 A - - sin^ 2^12 sin^ A 
- 2 sin^ 26'23 sin^ A 

+ 2 a Si3 sin 26'i2 sin 26'23 (sin 5cp sin A — cos 26'23 cos 5cp cos A) A sin A . 



(38) 



(39) 



4 Series expansion up to first order in a 

In this section, we expand the probabihties up to first order in the small parameter a while 
keeping their exact dependence on ^13. One expects these formulas to be useful for relatively 
large values of S13. In addition, they will correctly account for the "atmospheric" resonance 
driven by the parameters Amg^ and ^13. 

4.1 Matter of constant density 



a c 



12 ) 



The eigenvalues of the Hamiltonian (6) up to first order in a are 
El 

E, ~ ^a + A- C13) + ^as'i^ (Ci3 + 1 - ^cos2^i3) 

(Ci3-l + Acos2^i3) 



2E 
2E 



^ AmoT 

E3 ~ 

^ 2E 



where 



Ci3 = ^sin^ 2^13 + (>i- cos 2^13)' 



(40) 
(41) 

(42) 
(43) 



To calculate the probabilities to first order in a we followed two different approaches: First, 
we used the Cayley-Hamilton formalism as described in Appendix A. 2, to series expand the 
evolution matrix, and second, we considered the constant-density limit of the expansion of 
the evolution equation described in Appendix B.l. Both methods gave the same results, 
which is a useful cross-check of our calculations. Writing 



al3 



(44) 
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we 



obtain the following expressions for the Ve channel: 



p(o) = 1 

ee ^ q2 



sin' 2^13 . ^ 
— sm 613A, 



p(l) _ 



13 



X 



cosCisA sin Ci3 A cos 26'i3 - A 
— (1-^ cos 2^13) - A—- 



Similarly, for the i/, 



C'la 

channel we find 



d(o) 2 si'^^ 2^13 . 2 r< A 

^eM = ^23 —7^2 C13A , 

<^13 

p(i) o 2 2 ^in' 26' 

<^13 



2 2 sin2 2^i3 . 
' " — -7^ sm Ci3 A 



X 



COSC13A sin Ci3 A cos 26*13 - A 
— (1 cos 2^13) -^—^^ 

<-^13 <-^13 <-^13 



+ Si3 sin 26*12 sin 26*23 ^^".^2^^ | sin^cp [COSC13A - cos(l + ^)A] C13 

+ cos 5cp [Ci3 sin(l + A) A - (1 - A cos 26*13) sin C13A] | , 
and finally, for the — > v,,. channel we have 



^i? = ^ sin' 2^23 



1 - 



cos 



^pi)sm'i(l + X-C,3)A 



( ^ COS ^i/i 



Ism -(l + ^ + Ci3)A- -— ^2 — sm C13A 



P« = -^ sin' 2^23 a|2 



^2 
^13 

■' -' ^(1 + 25^3^ + ^') cosCi3Asin(l + A)A 



2 _ 2 2 ^ I 

^12 *12'5l3^2 ' 
•-^13 



22 22 22 /'22 2\/i1 ^^-^ ^^13^^ / /i \ a 

'12^13 ~ <^12'^13 + '^12*13 + i*12*13 ~ ^12 j 7^ C0S(1 + AjA 

<-^13 



-r ^ [1-12^-13 ^"12*13 

sin^ 26*13 sin C13 A 



X 



"12 ^2 

^13 ^13 

A. . ,^ . 74cos2^i3 — A . ^ . . „^ ^ ^ . 

— sm(l + A) A + — smCigA - (1 - A cos 2^13) cosCisA 

A A Ci3 

si3 sin 2^12 sin 2^23/^ 2 • c sinCi3A 
77T^ i 2ci3sm5cp — 7^ [cos(:;i3A - cos(l + A)A\ 

^^^13 Oi3 

- COS 26*23 COS 5cp(l + A) [cos Ci3 A - cos(l + A)Af 
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+ COS 2^23 COS 6cp 



• /-. ^\ A cos26'i3 — ^ . „ . 
sm(l + A) A H ^ sm C13A 



X 



C 



Ci3 

(1 + A) sm(l + A) A 



13 



(50) 



In the above formulas, one may identify the effective mixing angle in matter in the 1-3 sector 
^^3, which is determined through the expressions 



sin 2e[ 



sin 2^13 



13 



C 



cos 29', 



cos 2^13 - A 



13 



13 



c 



13 



(51) 



appearing frequently in Eqs. (45)-(50). Furthermore, the combination C13A appearing as the 
argument of sine or cosine corresponds to the effective A in matter. In the limit when ^13 is 
small, one has C13 ~ ^4 — 1, and expanding Eqs. (45)-(50) up to second order in S13 yields 
the double expansions given in Sec. 3, except for the terms of order a^. 

Equations (40)-(42) and (45)-(48) have previously been derived in Ref. [35] and confirmed 
in Ref. [36]. Expressions (49) and (50) are new. 

4.2 Vacuum oscillation probabilities up to first order in ct 

Taking the limit A ^ in Eqs. (45)-(50), it is straightforward to obtain the neutrino oscil- 
lation probabilities in vacuum to first order in a: 

pa)vac _ ^^^2^ g-j^2 2^^3 gjj^ 2A , (53) 

pf)^^^ = sin^ 2^135^3 sin^ A , (54) 

pWva. ^ _^^2^ 3i^2 2^^3^2_^ gin 2 A 



pWvac ^ ^4^ g-^2 2^23 Sin^ A , 



(55) 
(56) 



+ A sin 2^^i2Si3C^3 sin 26'23(2 sin (5cp sin^ A + cos (5cp sin 2 A) , 
-^13 sin" 26^23 r^"^ 

pWva. ^ _^^2^ 3i„2 2^^3(^2^ _ ^2^^2^) 2A 

-I- A sin 2^12 Si3C^3 sin 2^23 (2 sin bcv sin^ A — cos 5cp cos 2^23 sin 2 A) . (57) 



5 Series expansion up to first order in sis 



In this section, we expand the probabilities up to first order in the small parameter S13 while 
keeping their exact dependence on a. These formulas are expected to be useful whenever 
neutrino oscillations driven by the solar mass squared difference Am|]^ are important. In 
practice, this means low neutrino energies and long baselines. 
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5.1 Matter of constant density 



The eigenvalues of the Hamiltonian (6) up to first order in S13 are given by 

'A 



El ~ ^ 

2E 



a 



2E 



2 + 2(1-0,.) 



E3 ~ ^ 

^ 2E 



where 



Ci2 = \ sin^ 2^12 + cos 2^12 - 



A 



a 



(58) 

(59) 
(60) 

(61) 



Note that these eigenvalues are independent of the expansion parameter S13, which is consis- 
tent with Eqs. (27)-(29), where the lowest-order in S13 corrections to the eigenvalues appear 
only at order s^g. 

As in the case of the single expansion in a, we have calculated the probabilities using 
two different methods: by series expanding the evolution matrix using the Cayley-Hamilton 
formalism as described in Appendix A. 2, and using the constant-density limit of the single 
expansion in S13 of the evolution equation described in Appendix B.2. We have found identical 
results. Writing the probabilities as 



Pr 



we obtain the following expressions for the i/e i^e channel: 

sin^ 2^12 



p(0) ^ I _ 
ee 

Pii) = . 



fi2 



sin Q!Ci2A , 



(62) 

(63) 
(64) 



The absence of any first order corrections to Pge is consistent with Eq. (30). Similarly, for 



the Up 



Vfj^ channel we find 



siir 



20 



12 



'23 



sin^ aCi2 A , 



(65) 



12 



p(i) 



sin 2^12 . 



Cv2 

— cos 5cp 



sin 2^23 



(1 — a) sin q;Ci2A 
I- A-a + Aacl 



12 



sin 5cp [cos q;Ci2 A — cos(74 -|- o; — 2) A] 



, cos 2^12 - ^ 
sin(A + a - 2) A - sin q;Ci2 A ' 



aACi2 sm^2e 



12 



C 



12 



2(1 -a) 



C2 



12 



(66) 
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and for the — > Vr channel we have 



P(? = -Sin2 2^. 



23 



1 - - 



1 sin^ 29 



12 . 2 



sin^ aCuA - cos(q;Ci2 + A + a- 2) A 



'12 



COS 2^12 - - , 

- I 1 ^ sin aCi2A sm{A + a-2)A 



C 



12 



(67) 



sin 2^12 . 
-PAZ'' = — sin26'23 



X 



Ci2 

aACi2 



I — A — a + Aac\2 



cos 2^23 COS 5cp 



(cos q;Ci2 A - cos(A + a - 2) A)' 



COS 2^ — 

+ I ^ ^ sin aCiaA + sin(A + a - 2) A 



X 



Cl2 

'cos2^^i2-^ , 2(1 -a) 



C 



+ 



12 



sin q;Ci2A + sin(A + a - 2) A 



12 



+ sin 5cp (1 — ck) (cos q;Ci2A — cos(74 + a — 2) A) sin q;Ci2A > . 



(68) 



In this case, one may identify in Eqs. (63)-(68) the effective "solar" mixing angle in matter 
6*^2, which is determined by 



sin2^; 



sin 26' 



12 



12 



c 



12 



COS2^12-^ 
COS 2^12 = 77^^ 



(69) 



'12 



The combination aCi2A appearing as argument of sine or cosine corresponds to oscillations 
with the "solar" frequency in matter. Furthermore, we note that in the limit when a is small, 
one has C12 ~ A/o; — cos 26*12, and expanding Eqs. (63)-(68) up to second order in a yields 
the double expansions given in Sec. 3, except for the terms of order sfg. 

Neutrino oscillations probabilities in matter of constant density expanded to first order 
in Si3, but exact in ct, presented in this subsection, have not been previously published and 
are entirely new. 

5.2 Vacuum oscillation probabilities up to first order in S13 

Taking the hmit A — > in Eqs. (63)-(68), one recovers the vacuum probabihties expanded 
to first order in S13: 



pWvac _ ^ _ gi^2 2^^^ sin2 ^ 
pa)vac ^ ^ 

pWvac ^ sin2 2e^2cl^ sin^ aA , 

pW^ = cos 5cp sin 2^12 sin 2^23 [sin^ A - sin2(l - a) A + cos 2^12 sin=^ a A] 



(70) 
(71) 
(72) 
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+ - sin Sep sin 2012 sin 26'23[- sin 2A + sin 2(1 - a) A + sin 2q;A] , (73) 

pWvac _ gi^2 2^^^ gij^2 ^ ^ ^2^ gj^2 2e^^\^m^{\ - a) A - sin^ aA] , (74) 
pWvac ^ g-^ 26*12 sin 26*23 cos 26*23 [- sin^ A + sin2(l - a) A + cos 2^i2 sin^ a A] 

+ \- sin 5cp sin 26^12 sin 26'23[- sin 2A + sin 2(1 - a) A + sin 2aA] . (75) 

6 Qualitative discussion and tests of accuracy 

In this section, we discuss the qualitative behavior of the neutrino oscillation probabilities 
and the relevance of matter effects (Sec. 6.1). We also assess quantitatively the accuracy 
of the various expansions and compare each type of expansion with an exact numerical cal- 
culation of the corresponding probability (Sees. 6.2 and 6.3). We examine in detail which 
type of expansion is most accurate, depending on the values of the fundamental parameters 
(especially a and S13), and on the experimental configuration characterized by the neutrino 
energy E and the basehne length L. In Sec. 6.4, we discuss issues related to the applica- 
tion of our formulas to neutrino oscillation experiments. As we have shown in Sec. 2, the 
probabilities for all oscillation channels can be obtained from the expressions for Pg/^ and P^^ 
[sec Eqs. (16)-(22)]. Therefore, we focus in the following on the Vf, z/^ channel; we briefly 
comment also on the other neutrino oscillation channels. 

6.1 The relevance of matter effects and the probabihty Pep, 

Before considering the accuracy of the expansion formulas, we discuss in this subsection 
the features and the relevance of matter effects for the Vf. — > i/^ oscillation probability Pg^^. 
The discussion applies also to Per; matter effects on the probabihties P^^, Ptt, and P^^,- are 
significantly weaker than those on Pg^ and Pg^. Figure 1 shows contours of Pg^ for matter 
of constant density calculated numerically from Eq. (6) without any approximations, for a 
wide range of baseline lengths and neutrino energies. Many features of this figure can be 
understood by considering the expression for the two-flavor neutrino oscillation probability 
in matter of constant density: 

^•2 c\r\ j 

P = sin^ CA2. , C = ^sin^ 29 + (cos 29 - A^.f , (76) 

where, in analogy to Eqs. (25) and (26), we define = Am^L/{AE), A2U = 2EV/ Arn^, 
and 9 and Am^ are the generic two-flavor neutrino oscillation parameters. From Eq. (76) 
it is clear that for A21, <^ 1 one has C ~ 1, and vacuum neutrino oscillations with P ~ 
sin^ 2^ sin^ A2,y are recovered. For A2V — cos 2^ one can see in Eq. (76) also the MSW 
resonance [25], which leads to C = sin 2^ in Eq. (76) and to the effective mixing angle in 
matter sin^ 29' = sin^ 29 /C^ = 1. The general resonance conditions for three flavors are much 
more complicated than those in the two-flavor case. However, due to the hierarchy of the 
mass squared differences and smallness of 6'i3, one can use, to a very good approximation, the 
respective effective two-flavor picture. This leads to the "solar" and "atmospheric" resonance 
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IZZ1<0.1% 0.1% -1% P"«l%-5% ™5%-10% ™>10% 




E [GeV] 

Figure 1: Contours of Pg^t calculated numerically without approximations for sin^ 2^13 = 
0.02, Am|i = 2 ■ 10-3 eV^ a = 0.03, 612 = 33°, 623 = 45°, and 6cp = 0. Matter of constant 
density p = 3g/cm3 is assumed, and the probabihty is averaged over a Gaussian energy 
resolution of 1%. 



energies shown as vertical lines in Fig. 1 



A 2EV . 2EV 



The "atmospheric" resonance is clearly visible in Fig. 1 at Ey-^s — 8.6 GeV and L > 5500 km. 
The "solar" resonance occurs at the energy Ej-^s — 0.11 GeV. Note that the solar resonance is 
not as pronounced in Fig. 1 as the atmospheric one, since the neutrino oscillation amplitude 
in vacuum sin^ 26*12 is already quite large. 

Far enough to the left of the vertical lines in Fig. 1 marking the solar and atmospheric 
resonance energies one has vacuum neutrino oscillations with the "solar" parameters 612, 
Am2i and with the "atmospheric" parameters 613, Amg^, respectively. Indeed, the typical 
L/E pattern is clearly visible in Fig. 1 to the left of the vertical lines. The diagonal lines in- 
dicate the constant values of L/E corresponding to the first solar and atmospheric oscillation 
maxima in vacuum, given by the conditions 

aA = ^ = - , A = ^ = - , 78 

iE 2' AE 2' ^ ' 

respectively. To the right of the vertical lines the probability is dominated by matter effects. 
For A2U ^ 1 we have C ~ A2V in Eq. (76), which has two consequences: First, since 
A2vA.2v = VL/2, the oscillation frequency becomes independent of energy, and second, the 
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oscillation amplitude becomes suppressed at high energies because sin^ 26/ A^^^ oc 1/ E"^. Both 
these effects are apparent in Fig. 1 for neutrino oscillations with the solar as well as the 
atmospheric frequency. 

In addition to the dependence on the neutrino energy, matter effects depend crucially on 
the baseline. Prom Eq. (76) one can see that the two-flavor oscillation probabihties approach 
the vacuum ones for 0^2^ -C 7r/2. Far above the MSW resonance energy, i.e. for 3> 1, 
one finds C ~ A2U. Hence, the short-baseline limit is C^2v — A.2u^2u = VL/2 <^ 7r/2, leading 
to vanishing matter effects for L ^ i^mat = tt/V ~ 5453 km. The horizontal lines in Fig. 1 
indicate this baselines of the first "matter effect maximum" for large energies at L = Lmat) 
and the first "matter effect minimum", L — 2L^at — 10907 km. At the baseline 2Li„at; where 
AA — VL/2 = TT, the terms proportional to and asi^ in the double expansion Eq. (31) 
disappear, and only the term proportional to 5^3 survives. Therefore, this baseline is especially 
useful to measure S13, since ambiguities due to parameter degeneracies are avoided [38]. A 
dedicated analysis of this "magic baseline" can be found in Ref. [46].^ Note that the distance 
-^mat depends only on the matter potential V, i.e. the energy E as well as neutrino masses 
and mixings do not enter. 

Let us now consider the energy region well below the MSW resonance energy, where 
A2U ^ 1- The short-baseline behavior of the oscillation probability in this region can be 
understood in the two-flavor picture by expanding Eq. (76) for flxed A21/ {i.e., flxed L/E) 
assuming A2U ^ 1 and A2uA2u <^ 1- One flnds 

pmatter _ pvacuum _ 3 sin^ 29 COS 29 siu Ag, (siu Ag, - A2, COS A2,) ^2. ■ (79) 

This relation shows that p™3,tter _ pvax;uum (decreases linearly with A21, oc E, which means 
that the matter effects vanish for very small energies, as mentioned above. In addition, for 
a flxed energy the right-hand side of Eq. (79) becomes zero for A21, -C 1 or, equivalently, 
L <S Lvac = A^T^E/ Arn?. Therefore, the two conditions for the smallness of matter effects due 
to the short baseline can be written for all energies as L ^ Lmin = min(Lvac, -t/mat)- 

The relevance of matter effects is illustrated in Fig. 2. In the upper left panel we show 
the probability difference \P^^^^^ — PJ^^^^^\ for the case a = 0, such that the oscillations 
due to the "solar" frequency are switched off and only matter effects are present which are 
related to the "atmospheric" resonance energy at 8.6 GeV. It can be seen that the difference, 
as well as the relative difference (lower left panel), vanish for L <^ I>min- The left plots of 
Fig. 2 show that the differences vanish also for small energies and long baselines. This can 
be understood as follows. For E -C ^^res the oscillation amplitude in matter is always close 
to that in vacuum. For very short baselines, this is also true for the oscillation phases. At 
intermediate baselines, the matter-induced correction to the oscillation phase, though much 
smaller than the main "vacuum" term, can become comparable to unity and so cannot be 
ignored. However, at long baselines L ^ Lyac the averaging regime sets in, which implies that 
any shifts in the oscillation phase become unimportant and the vacuum oscillations limit is 
regained. 

■^It should be noted that the quoted number for Lmat holds for a constant matter density of p = Sg/cm"^. 
It will differ for larger matter densities like those in the mantle of the Earth, or even more drastically 
in astrophysical or cosmological applications of neutrino oscillations. In addition, the approximation of 
constant matter density is then frequently not justified. For a realistic Earth matter density profile one finds 
2Xmat ^ 7250 km [46]. 
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Figure 2: Contours of IP^T**'' " ^"""1 (upper panels), and IP^T"'' - P^'^MPlp. 
(lower panels), where p™'^**'^'' (Pg™'^""™) is the exact oscillation probability in matter (vacuum) 
for sin^26'i3 = 0.05, = 2 ■ lO^'^eV^, and 623 = 45°. In the left panels a = 0, whereas in 

the right panels a = 0.03, 612 = 33°, and 6cp = 0. Matter of constant density p = 3g/cm^ 
is assumed, and the probabilities are averaged over a Gaussian energy resolution of 2%. The 
straight lines indicate the relevant resonance energies, first oscillation maxima, and the first 
matter effect maximum (see also Fig. 1). 



The above discussion describes the regions where | p™a,tter _ pvacuum | vanishes or becomes 
small for a = in the left plots of Fig. 2. One can also understand in this way why the 
difference \P^'^^^'^^ — PJ^™*^™! disappears more slowly towards low E and L along the lines 
of constant L/E in the upper left plot of Fig. 2. The structures seen in the plot emerge 
from the matter dependent phase shifts of oscillation probabilities with more or less equal 
amplitudes. This behavior can be understood in the two- flavor picture from Eq. (79), which 
implies that p^^tt^r _ pvacuum gj-Q-^g linearly with E along the lines of constant A21, or, 
equivalently, constant L/E. This growth is modulated in the A2i/-direction by the factor 
sin (sin A2jy — A2;y cos A2;y), which describes nicely the details of the structures extending 
along lines of constant A21, towards lower energies. 

The right plots of Fig. 2 show the probability difference \P^^^^^^ — P^/^'^'^"™! and the 
relative difference for the case a = 0.03. This figure contains the structures stemming from 
the atmospheric resonance energy, which were already present in the a = case (left plots). 
In addition, similar structures show up around the solar resonance energy, reaching to lower 
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Figure 3: Rows 1, 2, and 3: absolute errors of the three types of expansions for P^,^ as 
functions of neutrino energy E and basehne length L, for Am|]^ = 2 ■ 10~^eV^, a = 0.03, 
6^12 = 33°, 6^23 = 45°, 6cp = 0, and several values of sin^26'i3. Matter of constant density 
p = 3g/cm^ is assumed. Row 4 shows contours of P^^ calculated numerically. 

energies and, due to the larger effective mixing angle, also to shorter baselines. The right 
plots of Fig. 2 show that in the realistic case of two mass squared differences matter effects 
are rather important and have to be taken into account in a large domain of the physically 
interesting parameter space. 

6.2 Comparing the accuracy of the three types of expansions 

In order to test the accuracy of our analytic expressions, we shall now compare the values 
of Pe/i obtained from the expansion formulas Eq. (31) for the double expansion, Eqs. (47) 
and (48) for the single expansion in a, and Eqs. (65) and (66) for the single expansion in 
Si3 with the exact values of Pe^ calculated numerically. Figure 3 shows the absolute errors 
ipexpansion _ pexact | ^ ^^^^ pig. 4, the relative errors |pe^pansion _ pexact|/pexact ^^le three 

types of expansions as functions of the neutrino energy E and baseline length L for various 
values of sin^2^i3. For reference, we display in the lowest row of graphs in each of these 
figures also the probability Pe^ itself, calculated numerically without any approximations. 

First, we observe that in general the absolute errors shown in Fig. 3 are rather small — 
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Figure 4: Rows 1, 2, and 3: relative errors of tlie tliree types of expansions for P^.^ as 
functions of neutrino energy E and baseline length L, for Am|]^ = 2 ■ 10~^eV^, a = 0.03, 
6^12 = 33°, 6^23 = 45°, 6cp = 0, and several values of sin^26'i3. Matter of constant density 
p = 3g/cm^ is assumed. Row 4 shows contours of P^^ calculated numerically. 

at the 0.1% level in a large part of the parameter space, whereas the relative errors shown 
in Fig. 4 are considerably larger, due to the smallness of the probability itself. Next, we 
note that the series expansions in a, i.e., the double expansion in a and sis and the single 
expansion in a, are only valid for 

aA = ^^^^ < 1 , or ^ < 10^ km/GeV , (80) 

i.e., far below the first solar maximum. The obvious reason is that expanding terms of the 
type sinaA is only valid if aA is small, and hence, these two types of expansion cannot 
account for neutrino oscillations with the "solar" frequency. Note that for S13 = the single 
expansion in a gives Pe/x = 0? since in that case the lowest-order term of P^^ is proportional 
to a^, and our single expansion in a only contains terms up to first order in a. This explains 
why at very small values of S13 the double expansion (which includes terms of second order 
in a) is more accurate than the single expansion in a. In contrast, the neutrino oscillations 
with the "solar" frequency (for which the condition in Eq. (80) is violated) are in general 
very well accounted for by the single expansion in S13, since it retains the exact dependence 



20 



of the probability on a. This expansion is the best one for relatively small values of S13 and 
large values of L/E. 

Figures 3 and 4 allow us to put the above observations on a more quantitative basis. In 
the region oi L/E given in Eq. (80) (lower-right parts of the graphs), where the oscillations 
are mainly driven by the "atmospheric" mass squared difference Amg^ , the double expansion 
and the single expansion in a work rather well. For not too large values of sin^ 2^13 the 
double expansion is most accurate, with absolute errors at the 0.1% level and relative errors 
not exceeding 1% (5%) for sin^ 2^^i3 = 0.001 (0.01). However, for values of sin^ 2^^i3 close to 
the current upper bound the single expansion in a becomes better, with relative errors smaller 
than 1%. The reason for this is that for large values of sin^ 2^13 the neutrino oscillations driven 
by the "atmospheric" frequency and mixing angle ^13 completely dominate the probability, 
and, in addition, for relatively large values of L and energies close to 10 GeV, the atmospheric 
resonance becomes important. Since the single expansion in a is exact in S13, it describes the 
case of relatively large S13 well, and the atmospheric resonance is also correctly accounted 
for. At the same time, the accuracy of the double expansion becomes slightly worse near 
the first atmospheric maximum and the resonance. As was pointed out in Ref. [36], for large 
values of S13 the accuracy of the double expansion can be improved by replacing the term 
proportional to sf^ in Eq. (31) by the term given in Eq. (47) as the zeroth order term in the 
single expansion in a. For a = this term describes the probability Pen exactly to all orders 
in S13. 

As can be seen from Figs. 3 and 4, the single expansion in S13 gives a rather poor de- 
scription of the region of L/E defined in Eq. (80). The reason is that with only terms of first 
order in S13 it is not possible to obtain a correct description of neutrino oscillations driven 
by Am|^ and ^^13. This is also reflected by the fact that the lowest-order in S13 terms in 
the eigenvalues of the Hamiltonian are C(si3) [see Eqs. (27)-(29)]. On the other hand, the 
single expansion in 613 is rather accurate for large values oi L/E violating the condition (80) 
(upper-left parts of the graphs) and relatively small values of ^13: for sin^ 2^13 < 0.001 the 
accuracy is typically better than 1%. 

To conclude this subsection, we show in Fig. 5 which type of expansion provides the most 
accurate expression for Pe^, depending on the values of the expansion parameters a and S13 
and for a number of fixed values of neutrino energy E and baseline length L. These plots 
change very little when the fundamental neutrino parameters Am§j^, 6*12, 6*23, and 6cp are 
varied within their allowed ranges, and also when one switches over to the other neutrino 
oscillation channels. As expected, one observes as a general trend that the single expansion 
in a is best for small a and large S13, whereas the single expansion in S13 is best for small 
Si3 and large a. The double expansion is most accurate in a region where a and S13 are of 
comparable order of magnitude. 

In agreement with the discussion related to Figs. 3 and 4, we find from Fig. 5 that in the 
low-energy regime ~ 0.1 GeV the single expansion in S13 is most accurate in almost the 
entire a-9i3 plane. For values oi L/E satisfying the condition (80) and energies larger than a 
few GeV the double expansion is most accurate in a large fraction of the a-Oi^ plane; only for 
values of S13 close to the current upper bound does the single expansion in a become better. 
Furthermore, we note that for E ^ GeV and L > 5500 km the atmospheric resonance is 
important, which leads to a better accuracy of the single expansion in a (see the rightmost 
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Figure 5: The plot shows which type of expansion for Pe^ is most accurate, depending on 
a and S13 and for different values of neutrino energy E and baseline length L. The values 
of the fundamental neutrino parameters used are Am^i = 2- lO-^eV^, 012 = 33°, 023 = 45°, 
and 6cp = 0. Matter of constant density p = 3g/cm^ is assumed. The dotted lines indicate 
the 3a allowed range of the parameter a. 



panel in the third row of Fig. 5). 

We conclude that the double expansion in both a and S13 is most accurate in a wide 
region of the parameter space, where a and S13 are roughly of the same order of magnitude. 
The single expansion in S13 has to be used whenever neutrino oscillations with the solar 
frequency are important {e.g., in the low-energy regime), whereas the single expansion in a 
is most accurate for values of S13 close to the upper bound, or in cases where the atmospheric 
resonance is important. 

6.3 Accuracy of the double expansion 

Motivated by the fact that the double expansion is most accurate in a wider region of the 
parameter space than the single expansions, we present in this subsection some more accuracy 
tests for it. Figures 6 and 7 show the relative errors for Pg^ given in Eq. (31) and P^r given in 
Eq. (34), respectively, for neutrinos, antineutrinos, and normal and inverted hierarchies. We 
note that the relative accuracy of Per and 1 — P^e is very similar to the one of Pg^ (Fig. 6), 
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E [GeV] E [GeV] E [GeV] E [GeV] 



Figure 6: Relative errors of the double expansion for for neutrinos (rows 1 and 2), and 
Pep for antineutrinos (rows 3 and 4), for normal hierarchy (NH) and inverted hierarchy (IH). 
The values of the fundamental neutrino parameters used are [Amg^l = 2-10-3eV^ a = 0.03, 
6*12 = 33°, 6*23 = 45°, 6cp = 0, and several values of sin^26'i3. Matter of constant density 
p = 3 g/ cm^ is assumed. 

and the accuracy of P^^ and P^-t is similar to that of P^^ (Fig. 7). 

In the region of L/E defined in Eq. (80), the accuracy of Pe^ is roughly between 1% and 
5%. It should be noted that the absolute errors of the double expansion of Pe^ depend very 
weakly on the choice of the fundamental neutrino parameters within their allowed ranges, and 
they are always very small, similar to the first row of plots in Fig. 3. However, the relative 
errors shown in Fig. 6 are rather sensitive to the fundamental parameters values in the region 
of low L/E. The reason is that the probability P^^ itself is typically very small in this region, 
being as tiny as 10~^, or even smaller. Hence, the relative error, being a ratio of two small 
numbers, is quite sensitive to variations of the parameters (see, e.g., the second column in 
Fig. 6). In contrast, P^^ is rather large, because of the oscillations in "atmospheric" channel 
driven by the large Amg^ and maximal or nearly maximal mixing sin^26'23 — 1. This is also 
clear from Eq. (34), since P^^ has a term of zeroth order in both a and S13, while Pg^ from 
Eq. (31) is of second order in these parameters. Hence, in the case of P^r, the relative errors 
are of the same order as the absolute errors, which are always very small. The errors shown 
in Fig. 7 are smaller than 0.1% for sin 2^13 < 0.01 and smaller than 1% for sin^ 2^13 < 0.1. 
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Figure 7: Relative errors of the double expansion for P^^ for neutrinos (rows 1 and 2), and 
PjPf for antineutrinos (rows 3 and 4), for normal hierarchy (NH) and inverted hierarchy (IH). 
The values of the fundamental neutrino parameters used are [Amg^l = 2-10-3eV^ a = 0.03, 
6*12 = 33°, 6*23 = 45°, 6cp = 0, and several values of sin^26'i3. Matter of constant density 
p = 3 g/ cm^ is assumed. 

6.4 Probability expansions and relation to experiments 

The discussion of the accuracy of the different expansions can be used to identify the best 
set of equations for a given neutrino oscillation experiment. The first question is, however, 
whether matter effects are relevant or if the much simpler vacuum probabilities can be used. 
The absolute and relative differences between matter and vacuum probabilities vanish, as 
discussed in Sec. 6.1, for sufficiently low energies or for short enough baselines. The discussion 
showed that in experimentally interesting cases matter effects are almost always relevant for 
the ranges of L and E that we considered, li L/E is such that oscillations can be observed, 
then for fixed L/E the matter effects are negligible only at very small E (and consequently 
small L). An example where matter effects can be ignored at the percent level is given by 
the reactor disappearance experiments with energies of a few MeV and baselines up to a 
few kilometers. 

In order to identify the best expansion for a given experiment, one must take into account 
that the analyses of real experiments are based on event rates and not on probabilities. It 
is useful to distinguish here between "point sources" (accelerator beams, reactors, the Sun, 
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supernovae, etc.) where the observed neutrinos come from one point-hke region in space, and 
"extended sources" {e.g., the atmosphere) where the observed neutrinos come from many 
different positions in space. An event rate based analysis of point sources imphes essentially, 
for a given energy, an additional factor in order to account for the reduction of the 

flux as a function of the baseline L. In addition, the detection cross section is typically 
proportional to E^, where 1 ^ r < 2, depending on the energy range and on the detection 
process. In the two-flavor approximation, when going from the first oscillation maximum at 
^2u = 7r/2 to the nth oscillation maximum at A21, = {n — 1/2) tt this implies that the event 
rates drop by a factor l/(2n — 1)^ if the basehne is increased, or by a factor l/(2n — 1)'' if the 
energy is decreased. Hence, e.g., for a fixed energy the rate in the second (third) maximum 
is already reduced by a factor 1/9 (1/25). This explains why most current proposals for long 
baseline oscillation experiments aim at the first oscillation maximum of the "atmospheric" 
oscillations. The oscillations driven by the "solar" frequency are then a sub-leading effect, 
and in most cases the double expansion works very well. The overall relative precision of the 
double expansion is then typically in the range of a percent, or at worst, up to a few percent 
in this case. However, if 613 is close to its current upper limit, the single expansion in a is 
even more precise than the double expansion. On the other hand, if ^13 is tiny, well below 
the sensitivity limits of any planned accelerator experiment, the single expansion in S13 is 
more precise than the double expansion, especially for low energies. Note, however, that the 
double expansion has often a sufficient precision even when it is not the best expansion. 

The discussion of point sources does not change much for experiments which aim at 
higher oscillation maxima. The l/(2n — 1)^ fiux factor due to the beam divergence and/or 
the l/(2n — 1)*" factor due to the energy dependence of the detection cross sections reduce 
the observed event rates and limit all proposals to the first few oscillation maxima at best. 
The double expansion still works quite well for such proposals, as long as the sub-leading 
oscillations governed by the "solar" frequency stay in the linear regime, where an expansion 
in a makes sense. This leads to the condition 2n — 1 ^ ^ 40, which is fulfilled for the first 
few maxima. Higher values of n are currently not proposed, since the fiux and/or cross section 
would drop by a large factor, which would make the detectors and sources unaffordable. For 
largest (smallest) 613 the single expansion in a (sis) works again numerically even better then 
the double expansion, just like in the case of the first oscillation maximum. 

The discussion is somewhat different for sources which are not point-like. In that case, 
there is in general no 1 /L^ suppression of the flux, and the observed event rates are not dom- 
inated by the first oscillation maximum. An important example are atmospheric neutrinos, 
for which a wide range of baselines and energies contributes: the energy window of practical 
interest ranges from a 100 MeV to a few tens of GeV, while the basehne lies between 10 km 
and 10^ km. This wide range of energies and basehnes makes clear that none of the discussed 
expansions covers the full parameter range at the percent level. In general the double ex- 
pansion works quite well for small or moderate L/E, while the S13 expansion tends to work 
best for larger L/E values. At the same time, close to the atmospheric resonance or if the 
mixing angle ^13 is near its current upper bound and for not too large L/E, the a expansion 
is expected to be the most relevant one. 

In general, for all types of experiments and for any given values of L, E, and ^13 the most 
accurate expansion can readily be found with the help of Fig. 5. 
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7 Summary and conclusions 



In this paper we have presented three different sets of approximate formulas for the prob- 
abilities of neutrino oscillations in matter and in vacuum. We have shown that in general 
the probabilities for all possible oscillation channels for neutrinos as well as for antineutrinos 
can be obtained from just two independent probabilities by using unitarity and a symmetry 
related to the "atmospheric" mixing angle 62?, [see Eq. (14)]. One possible choice for these 
two probabilities is Pg^ and P^,-. We have derived the expressions for the neutrino oscilla- 
tion probabilities in matter of constant density, expanded in terms of the small parameters 
a = Am|^/Am3^, S13 or in both of them. Below we summarize the main features of these 
expansions: 

• Double expansion up to second order in a and S13 [Sec. 3.1, Eqs. (30) (34)] 

The neutrino oscillation probabilities are expanded in both small parameters up to 
second order. In general, these expressions are valid for aA = Am\^L/{AE) <^ 1 
or L/E <^ lO^km/GeV, i.e., if neutrino oscillations due to the solar mass squared 
difference Am|j are not important. The accuracy of the approximation is good for a 
wide range of the parameters. Typically, the relative errors of Pg^, Pgr, and 1 — P^e are 
between 1% and 5%, the relative errors of P^^^, Pj^r, and P^r and the absolute errors for 
all probabilities are of the order 0.1%. 

• Single expansion up to first order in ot [Sec. 4.1, Eqs. (45)-(50)] 

The neutrino oscillation probabilities are expanded in a, but the exact dependence on 
Si3 is retained. Like in the case of the double expansion, these formulas are valid in 
the region where the oscillations driven by the solar mass squared difference are not 
important: aA = Aml^L/{AE) < 1, or L/E < lO^km/GeV. The accuracy is better 
than the one of the double expansion for values of S13 > 0.1 close the current upper 
bound, or if the atmospheric resonance is important. For instance, for a matter density 
p = 3 g/cm^ and Am|i = 2 • 10"^ eV^ this is the case for £^ ~ 10 GeV and L > 5500 km. 

• Single expansion up to first order in S13 [Sec. 5.1, Eqs. (63)-(68)] 

The neutrino oscillation probabilities are expanded in S13, but the exact dependence 
on a — and therefore on Am\^ — is retained. Hence, these expressions are useful in the 
region where the oscillations due to the solar mass squared difference are relevant: 
aA = Am|j^L/(4£') > 1, or L/i? > 10*^km/GeV. In particular, this is the case for low 
energies ~ 0.1 GeV and L > 10^ km. For values ol L/E outside the "solar" regime 
this type of expansion is only useful for very small values of S13 < 10~^ 10~^. 

• Vacuum limit of the expansions [Sees. 3.2, 4.2, 5.2] 

The vacuum limit for each type of expansion is valid if, in addition to the requirements 
ensuring the validity of a given type of expansion, matter effects can be neglected (see 
discussion in Sec. 6.1). This usually imphes low energies or short basehnes. We find 

that in many realistic cases matter effects are of the order of a few percent. Hence, 
if an accuracy at that level is required, one should employ the formulas for neutrino 
oscillations in matter. 
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In conclusion, we have presented a collection of formulas for three-flavor neutrino oscilla- 
tion probabilities by deriving expansions in small parameters. We have performed a detailed 
analysis of the accuracy of these expansions and determined the parameter regions where 
they are most accurate. The expansions of the neutrino oscillation probabilities in matter 
of constant density are useful for the analytical understanding of the physics of future neu- 
trino oscillation experiments. Furthermore, we have also presented expansion formulas for 
the neutrino oscillation probabilities in arbitrary matter density profiles (see Appendix B), 
which can be apphed to a large class of problems. 
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A Hamiltonian diagonalization approach 

In this appendix, we present the details of the approximate diagonalization of the effective 
Hamiltonian of the neutrino system in the case of matter of constant density. 

A.l Diagonalizing the Hamiltonian perturbatively 

In order to derive the double expansions given in Sec. 3, we write the Hamiltonian of Eq. (6) 
as 

/\?T7 ^ 

Hc^^OMMUlOl, (Al) 

where M = O13O12 diag (0, a, 1) 0^2^13 + diag (^4, 0, 0), and the matrices and Us have 
been defined after Eq. (1). The matrix M can be explicitly written as M = (Amli/2E)~^H' 
by setting 6cp = in H', which is given in Eq. (Bl) below. First, we diagonalize the matrix 
M by M = WMW^ with M = diag(Ai, A2, A3) and W being a unitary diagonalizing matrix. 
This diagonalization is performed by using perturbation theory up to second order in the 
small parameters a and S13, i.e., we write M = M(°) + M^^) + M^'^\ where M^^) (M^^)) 
contains all terms of first (second) order in a and S13. One finds 

M(°) = diag(/l, 0, 1) = diag(AS°\ A^ , A^ )) , (A2) 
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asj^ «Sl2Cl2 Sl3 



M(^) = I asi2Ci2 OLc\^ I , (A3) 
Sis 



'13 



—OLS\2,s\2 



M(2) = I -asi3Si2Ci2 . (A4) 

-q;si3Si2 -asi3Si2Ci2 -s?3 / 

For the eigenvectors we write Vi = vf^ + v^^^ + v^'^^ . and since M'^'^^ is diagonal at zeroth order, 
we have vf'^ = Cj. Then, the first and second order corrections to the eigenvalues are given 

by 

^ = , (A5) 
and the corrections to the eigenvectors are calculated by 



(2) 



3+i ^3 



ej . (A8) 



The mixing matrix in matter is given by U' — O^^J^sW with W — (I'l, i'2, ^3), and the 
eigenvalues of the Hamiltonian [see Eqs. (27)-(29)] are obtained as Ei — [/^m\^/ {2E)\\ 
(i = l,2,3). 

In our paper, we do not in general order the eigenvalues according to their magnitude. 
Such an ordering would create problems as one would have to re-label the eigenvalues upon 
passing through each of the two MSW resonances. The ordering is actually unimportant if 
one is careful to assign the correct eigenvector to each eigenvalue. 

We also reiterate the point discussed at the end of Sec. 3.1: despite the fact that in the 
case of the double expansion the eigenvalues as well as certain entries of the leptonic mixing 
matrix in matter are divergent at A ^ and yl — > 1, the neutrino oscillation probabilities 
are finite in these limits. In particular, the correct vacuum probabilities are recovered in the 
limit ^4 ^ 0. The mentioned divergences arc of no concern to us, since we are interested 
in oscillation probabilities rather than in eigenvalues or the mapping between the mixing in 
matter and in vacuum, as was the case, e.g., in Ref. [36]. 
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A. 2 The Cayley-Hamilton formalism 



In order to derive a series expansion for the evolution matrix S in small parameters, we can 
use the following formula from Refs. [12, 14]: 

[E2E3l3-{E2 + E3)H + H^] 





{El 


— E2) {El 






Q-IE2L 




(E2 


— Ei){E2 


-Es) 




— Ei){E^ 


-E2) 



where H is the Hamiltonian in the flavor basis as given in Eq. (6) and Ei [i = 1,2,3) are 
the eigenvalues of this Hamiltonian. Inserting series expansions of the eigenvalues [see, e.g., 
Eqs. (40)^(42) and (58)-(60)] and the Hamiltonian H into Eq. (A9), it is straightforward 
to obtain an expression for S'^/j, and subsequently, for the neutrino oscillation probabilities 

-Pa/3 = |'S'/3aP- 



B Perturbative expansion of the evolution equation for arbitrary 
matter density profiles 

In this appendix, we present the details of the perturbative expansion of the neutrino evolution 
equation in the case of matter with an arbitrary density profile. We also give the formulas 
relevant for the particular case of matter of constant density. 

Consider the neutrino evolution matrix to) defined in Eq. (9). It is convenient to 
rotate it by the angle ^^23 according to Eq. (13). The evolution matrix in the rotated ba- 
sis S"(t, to) satisfies the evolution equation i(d/dt)5"(t, to) = H'{t)S'{t,to) with the initial 
condition S"(to,to) — 1 and the Hamiltonian H'{t) given by 



Sl2cls A21 + S^gAgi + V{t) S12C12C13A21 S13C13 (A31 - S?2A2l) e"''^'^^ 

H'{t) = I S12C12C13A21 CI2A21 -Sl2Ci2Si3e-'^^^A2l 

S13C13 (A31 - SI2A21) -Si2Ci2Sl3e'^^^ A21 C23A31 + S?2S?3^21 

Here 



(Bl) 



A21 = ^ , A3, ^ ^ . (B2) 
21 2E ' 2E ^ ' 

It is important to notice that H'{t) [and hence, S'{t, to)] does not depend on the mixing angle 
^23- This is a consequence of the specific parameterization of the leptonic mixing matrix U 
in Eq. (1), in which the matrix O23 is the leftmost one, and the fact that the matrix of the 
matter- induced potentials in the effective Hamiltonian H{t) commutes with 023- Once the 
matrix 5"(t, to) is found, one can use Eq. (13) to rotate back to the original fiavor basis. 
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It is convenient to decompose the Hamiltonian H'{t) as H'{t) = 11^(1) + Hj, where -f^o(^) 
is of zeroth order in some expansion parameter and Hj is the remainder. Accordingly, the 
evolution matrix can be written as S'{t,to) = SQ{t,to)S[{t,to), where 5*0(^,^0) satisfies 



i^5^(i, to) = K{t)S'^,{t, to) , Siito, to)^l. 



Then the matrix S[{t,to) satisfies 



i^^S[{t,to) = [S',{t,to)-'H'jS',{t,to)]S[{t,to) , S[{to,to) = 1. 



(B3) 



(B4) 



Up to now, no approximations have been made. Next, we shall find the evolution matrix 
perturbativcly. Let H[ be the part of Hj that is of the first order in the chosen expansion 
parameter. Then, to first order in this parameter, one finds 



S'{t,to)c^S'o{t,to)-iS'o{t,to) / [S'o{t',to)-'H[S',{t',to)]dt', 



(B5) 



to 



from which, upon rotation back to the original flavor basis, the evolution matrix S{t, to) is 
obtained. We shall now consider the appUcation of the above formahsm to the cases when 
the expansion parameter is either a or S13. 

B.l Expansion up to first order in a 



In this case, the zeroth-order Hamiltonian i^o(^) is obtained from Eq. (Bl) by taking the 
limit a — > {i.e., A21 0). The zeroth-order evolution matrix S'Q(t,to) can then be written 
as [26] 

/ u{t,to) v{t,to)e-'^^^ \ 

S'o{t,to)^ f{t,to) , (B6) 

\ -v*{t,to)d^^^ u*{t,to) J 

where u(t, to) and v{t, to) are to be found from the solution to the two-flavor problem governed 
by the mass squared difference Am|^, mixing angle 6*13, and matter potential V{t), and 



f{t,to) = exp 



i / A{t') dt' 

to 



Am|j^ 
2E 



+ V{t) 



(B7) 



The parameters u and v satisfy + \v\'^ — 1. The factor e ^^'^^ was pulled out from v for 
convenience (this will simplify some of the formulas given below) . 

The neutrino oscillation probabilities to first order in a in the case of matter with an 
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arbitrary density profile can be written as 

^ \u\'^ - 2 ^ s^2 l^[uv* (2 cos 2^13 7i + sin 2^13 h)] , 
aA 



(B8) 



aA 



Peix = S23l'"l + 2 S12 ^23 ^T^[uv* (2 COS 26*13 h + sin 26*13 -^2)] + sin 26*12 sin 26*23 

X {cos 5cY> lm[v*r{ci3h + S13/4)] + sin Scp Re[v*r{ci3h + S13/4)]} , (B9) 
P^^ = ] sin^ 2^23 [1 + l^r - 2Re(«/)] 



aA 



sin 2^23 \ Im 



+ 



- f*)v* ( cos 2^13 h + - sin 2^13 /2 



Ci2-f' - S12 ( Si3 -^5 + Ci3 h - ^ sin 26*13 h 



lm{uf) 



- ^ sin 26*12 sin26*23|cos5cp[c23lm[('u-/*)i;*(ci3/3 + Si3/4)] 

+ Im[(c23|M|^ - S23 - cos 26*23 Uf) (Ci3 I9 - Si3 Is)]] 

+ sin Scp [033 Re[(w - /*)v*(ci3 73 + S13 74)] 

- Re[(c23|M|^ - S23 - cos 26*23 uf) (C13/9 - S13/8) 
where the quantities Ii, I2, . . . , I9 are given by 

h = f u*{t',to)v{t',to)dt' , l2= I [v{t'Mf -u\t'Mf]dt' , 

U*{t',to)f{t',to)dt', 

f\u{t',to)\'dt', 

J to 



h = 



h = 



h = 



v{t',to)f{t',to) dt', 
f\v{t',to)\'dt', 

J to 



(BIO) 



Ir= [ [u(t', to)v(t', to) + u*{t', to)v*(t', to)] dt' , 

Jto 

h= fu*(t',to)nt',to)dt', h= fv(t',to)nt',to)dt'. 

Jto Jto 

Note that the integrals Jr,, Jg, and J7 are real, J5 + /g = t — to, and the integrals /g and Ig can 
be obtained from I3 and I4, respectively, by substituting / /*. 

In the case of matter of constant density, the parameters u, v, and / are given by 



, ^, Ci3Ai .cos26^i3-^ . CisAt 
u{t, U) = cos — h 1 sm 



C 



13 



.sin 2^13 . Ci3Ai 
v{t, 0) = -i—p^ sm — - — 

Ol3 Li 
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(Bll) 
(B12) 



/(i,0)=exp 



(B13) 



Inserting these expressions into the integrals Ii, I2, . . . , Ig and the obtained results into Eqs. (B8)- 
(BIO) leads to Eqs. (45)-(50). 

B.2 Expansion up to first order in S13 

This case was studied in detail in Ref. [28] (Appendix A), which we closely follow here. The 
zeroth-ordcr Hamiltonian HQ{t) is obtained from Eq. (Bl) by taking the limit ^13 — >■ 0. The 
zeroth-ordcr evolution matrix 6*0(^,^0) can be written as 



x{t,to) y{t,to) 

-50(^,^0) = I -y*{t,to) x*{t,to) 

g{t,to) 



(B14) 



where x{t, to) and y{t, tg) are to be found from the solution to the two-flavor problem governed 
by the mass squared difference Am^i, mixing angle 612, and matter potential V{t), and 



g{t, to) = exp 



-i/ A{t')dt' 



2E 



2E 



+ V{t) 



(B15) 



The parameters x and y satisfy + = 

The neutrino oscillation probabilities to first order in S13 in the case of matter with an 
arbitrary density profile can be written as 

|2 



\X\ 



3 lyp - sin 2^23 i^ygC) 



^23^23 \x- 9*\ + sin 26^23 Im[(a: - g*)g{s2sB - c^sD)] , 



where the quantities B,C, and D are defined as [47] 

B = B(t, to) = aly*,t(t, to) + hh^t{t, to) , 
C = C{t,to) = a*/:.,„(t,to) - h*Il,^{t,to) , 
D = D{t, to) = a*i;,^,^{t, to) + b*i:^,^{t, to) . 

Here the integral I^^s is defined as 

I^At,to)= [ ip{t',s)g{t',s)dt' , 

Jto 



and the quantities a and b are given by 



a — 
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2E 



Sl3(l-s?2«)e-^''«^ 



2E 



Sl3Sl2Cl2e' 



-Wcp 



(B16) 

(BIT) 
(BIS) 

(B19) 
(B20) 
(B21) 

(B22) 

(B23) 
(B24) 
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In the case of matter of constant density, the parameters y, and g are given by 

xit, 0) = cos + i ^"^'^^-^/" sin , (B25) 

.sin 2^12 . aCuAt 
y{t,0) = -i— sm , (B26) 

Ol2 Li 



g{t, 0) = exp 



.{A + a-2)A " 

1 ~ t 



(B27) 



Inserting these expressions into Eqs. (B19) - (B22) and the obtained results into Eqs. (B16)- 
(B18) leads to Eqs. (63)-(68). 
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